#include <cmath>
#include <cstdlib>
#include <vector>
#include <boost/math/special_functions/spherical_harmonic.hpp>

template<typename T>
class SHFunction
{
public:
	int nbands;
	std::vector<T> values;

public:
	static const T PI;

public:
	explicit SHFunction(int bands = 4);

public:
	static T random();
	static void random(T &theta, T &phi);
	static T realY(int l, int m, T theta, T phi);

	static void reverseAngle(T theta, T phi, T &rtheta, T &rphi)
	{
		rtheta = PI - theta;
		rphi = phi < PI ? phi + PI : phi - PI;
	}

	static void vectorToAngle(T x, T y, T z, T &theta, T &phi)
	{
		theta = std::acos(z);
		T cosphi = x / std::sqrt(x * x + y * y);
		phi = y >= T(0) ? std::acos(cosphi) : T(2) * SHFunction<T>::PI - std::acos(cosphi);
	}

	static void angleToVector(T theta, T phi, T &x, T &y, T &z)
	{
		T sint = std::sin(theta);
		x = cos(phi) * sint;
		y = sin(phi) * sint;
		z = std::cos(theta);
	}

public:
	T evaluate(T theta, T phi) const;
};
